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We review the statics and dynamics of electric-field domains on doped superlattices within a 
discrete drift model. A complete analysis of the construction and stability of stationary field profiles 
having two domains is carried out. As a consequence we determine the intervals of doping on 
which self-sustained current oscillations may appear under dc voltage bias. We have also studied 
the influence of doping, boundary condition and length of the superlattice on the self-sustained 
oscillations. Our study shows that there are bistability regions where either self-sustained current 
oscillations or steady states are reached depending on the initial condition. For a wide bias interval, 
the self-sustained oscillations are due to the formation, motion and recycling of electric-field domain 
walls inside the superlattice. There are biases (typically in the region of bistability) for which 
the strength of the high and low field domains changes periodically in time while the domain wall 
remains almost pinned on a few quantum wells. 
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I. INTRODUCTION 



Electrical transport in semiconductor superlattices 
(SL) has attracted much interest during the last years 
due to the very different interesting properties related 
to the artificial band structure. One of these features is 
the occurrence of stationary electric field domains which 
have been observed already in 1974 [1]. Due to advanced 
growing facilities and experimental techniques the com- 
plicated structure of the current-voltage characteristics 
which exhibits several branches, roughly equal to the 
number of quantum wells, could be resolved during the 
last decade [2-7]. In these experiments it was demon- 
strated that the different branches are connected to the 
formation of two domains with different values of the 
electric field inside the sample. Depending on the con- 
ditions, stable stationary field domains and traveling do- 
main boundaries may occur. In the latter case, the dy- 
namics of the electric field domains gives rise to time- 
dependent oscillations of the current. [8-11] 

During the last years it has been shown theoretically 
that the observed phenomena can be fairly well repro- 
duced by models which essentially combine the discrete 
Poisson-equation and rate equations for the carrier densi- 
ties in the different quantum wells [12-15]. Also the time- 
dependent current oscillations could be recovered in these 
models [14,16,17]. A prediction of spatio-temporal chaos 
in resonant-tunneling superlattices under dc+ac voltage 
bias has been made on the basis of the discrete drift 
model [18]. The influence of growth-related imperfec- 
tions on the SL behavior has been studied in [19]. These 
phenomena may also be described by discrete models 
with Monte Carlo dynamics incorporating single-electron 
tunneling effects (wich are important for slim superlat- 
tices and give rise to additional oscillations of the current 
[20,21]). 



In this paper we want to explain how these compli- 
cated phenomena are generated by such models. This 
provides a deeper insight into the basic mechanisms and 
helps to classify the results of various experiments and 
computer simulations. In particular we want to give an 
answer to the following questions: How is it possible to 
understand the appearance of the complicated structure 
of the current- voltage characteristic? What are the con- 
ditions for stability and oscillations and how can they be 
understood? What is the main difference to the Gunn 
diode, where hardly any stable domain states are ob- 
served? Where can the high-field domain be located with 
respect to the injecting contact? What are the mecha- 
nism^) explaining the current oscillations? How do dop- 
ing, electron velocity and boundary conditions influence 
the shape and frequency of the current oscillations? 

The paper is organized as follows: The model we use is 
described in the second section. The third section shows 
how the complex stationary current-voltage characteris- 
tic changes as the doping increases. In the fourth sec- 
tion we investigate the stability of the stationary states 
and prove an explicit criterion for the occurrence of sta- 
ble domain states. The fifth section provides a thorough 
numerical investigation of the oscillatory behavior gen- 
erated by the instability of the current discussed before. 
The last section contains our conclusions and the Ap- 
pendix is devoted to a proof that no self-sustained oscil- 
lations appear for low enough doping times number of SL 
periods. 



II. THE MODEL 

We consider a semiconductor superlattice where the 
lateral extension of wells and barriers is much larger than 
the total length of the SL, so that single-electron tun- 
neling effects (see e.g. Ref. [20]) are negligible. The 
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quantum wells (QW) are weakly coupled and the scat- 
tering times are much shorter than the tunneling time 
between adjacent QWs. Thus it makes sense to consider 
the electrons to be localized within the QWs and in lo- 
cal equilibrium at the lattice temperature. The current 
is mainly determined by the resonances between the dif- 
ferent energy levels in the QWs, which we denote by Ci, 
i = 1,2,..., in order of increasing energy counted from 
the bottom of the conduction band. For the biases of in- 
terest here, there are three important resonances C1C1, 
C1C2 and C1C3. If the intersubband relaxation is also 
fast with respect to the tunneling, in practice only the 
lowest subband is occupied. In this case it is a reason- 
able approximation to consider the QWs as entities char- 
acterized by average values of the electron density hi in 
the i-th QW and the electric field Ei between wells i and 
i+ 1, with i = 1, . . .N. The discrete drift model consists 
of the following system of equations for Ei , hi , and the 
total current density J(t) [16]: 



E i -E i - 1 = 2j(n i -N D ), 

dEi . ~ 

€— T + qv(Ei)ni = J, 
dt 

N 



(i) 

(2) 
(3) 



where i = 1, . . . , N . In this model Eq. (1) and Eq. (2) are, 
respectively, the one-dimensional Poisson equation (aver- 
aged over one SL period) and Ampere's law. e, No, I and 
q are the average permittivity, average donor concentra- 
tion, SL period, and the charge of the electron. v(E) is an 
effective electron velocity that has a peak at certain val- 
ues of the electric field connected to resonant tunneling 
C1C1, C1C2 and C1C3. In this paper we shall be con- 
cerned with phenomena occurring at fields higher than 
the first resonant peak C1C2, so that we shall omit the 
miniband peak C1C1 in our velocity curve; see the curve 
plotted in Fig. 1 by a dashed line. Phenomena at lower 
fields can easily be studied by adding the C1C1 peak to 
our curve, [23]. Equation (3) establishes that the average 
electric field is given by the dc voltage bias $. Notice that 
there are 2N+2 unknowns: Eo, E\, . . . , En, hi, . . . , njy , J 
and 2N+1 equations so that we need to specify one 
boundary condition for Eo plus an appropriate initial 
profile Ei(0). The boundary condition for Eo (the aver- 
age electric field before the SL) can be fixed by specifying 
the electron density at the first site, hi, according to (1). 
In typical experiments the region before the SL has an ex- 
cess of electrons due to a stronger n-doping there than in 
the SL [8-11]. Thus it is plausible assuming that there 
is an excess number of electrons at the first SL period 
measured by a dimensionless parameter c > — 1: 



hi = (1 + c)N D 



(4) 



c has to be quite small because it is known that a steady 
uniform-electric-field profile is observed at low laser il- 



lumination in undoped SL [8,14,16]. This observation 
allows us to infer the electron velocity directly from mea- 
sured current- voltage data [14]. Another possibility is to 
derive the electron velocity from simple one-dimensional 
quantum-mechanical calculations of resonant tunneling, 
as was done by Prengel et al [13]. They used a more 
complicated discrete model with two electron densities 
corresponding to the populations of the two lower en- 
ergy levels of each QW. Their model reduces to a form of 
ours when the large separation between the time scales of 
phonon scattering, resonant tunneling and dielectric re- 
laxation is taken into account. An earlier discrete model 
is due to Laikhtman [24] and an attempt at deriving a 
discrete model from quantum kinetics can be found in 
[12,15] (only miniband conduction was considered, not 
resonant tunneling). 

For the calculations that follow, it is convenient to 
render the equations (l)-(4) dimensionless by adopting 
as the units of electric field and velocity the values at 
the C1C2 peak of the velocity curve, v(E), Em and vm 
(about 10 5 V/cm and 427 cm/s, respectively, for the sam- 
ple of Ref. [10]). We set [16]: 



Ei 



Ek_ 

Em 
v 

vm ' 



_ qlhj 
cEm 

gN D v M t 
cEm 



J 



qN D v M 



N Em I 



(5) 



By using Poisson's equation to eliminate the electron 
density in favor of the electric field we obtain the di- 
mensionless equations: 

«1 = m - (i + ^r^) ( 6) 

for i=l,...N, 



1 N 



N 



(7) 



and the boundary condition 

E {t) = Ei{t)-cv. (8) 
The dimensionless parameter v, is defined by 

N D ql 



e Em 



(9) 



which yields about 0.1 for the SL used in the experiments 
[10,11]). The constant voltage condition (7) determines 
the current to be 



1 N ( 



Ei — Ei-i 



v(Ei) (10) 
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With the choice (5), the dimensionless velocity v(E) has 
a maximum at E = 1 with i>(l) = 1. Throughout this pa- 
per we use the function v(E) which is plotted in Fig. 1 by 
a dashed line. Besides having a maximum i>(l) = 1, it has 
a minimum at E m ps 1.667 with v(E m ) = v m ps 0.323. 
Nevertheless nearly all of the features discussed in the 
following are independent of the exact shape of the func- 
tion v(E). We only impose the restrictions that v(E) > 
for E > and the existence of a minimum at E m > 1 . 

For any <f> > 0, we shall assume that the initial electric 
field profile is strictly positive and that the electron den- 
sity is non-negative: Ei(0) > 0, n;(0) = E i (0)-E i - 1 (0) + 
v > 0, Vi. This is reasonable unless <f> is close to zero (but 
then the C1C1 peak in the velocity curve should be re- 
stored) or the boundary conditions are unrealistic. From 
the equations and our assumption on the initial field pro- 
file, it follows that E{ > and n 8 - = E{ — Ei-i + v > 
for all positive times. The model equations have interest- 
ing properties concerning the monotonic behavior of the 
electric fields with respect to the QW number. These are 
summarized in the following Lemmas: 

Lemma 1 If the fields of two adjacent QWs are identi- 
cal, i.e., Ek = Ek-i holds for some k with 2 < k < N, 
there is at least one i from 1 < i < k with Ei = Ei-i and 
d(Ei — Ei-\)/dt ^ 0. For c = there is the additional 
possibility that Eo = E\ = . . . = E k holds. 

Lemma 2 If c > and the field distribution is mono- 
tone increasing at t = (Ei(0) > £ , 8 '_i(0), Mi), it will 
continue being monotone increasing with respect to the 
well index i for all later times t > 0. (For c < 0, the 
same holds for a monotone decreasing field distribution). 

Both lemmas can be easily proved by using the prop- 
erty 

j t (E k - E k .,) = Ek ~^ k - X K^-i), (11) 

which holds for the particular k considered in the first 
Lemma, and the boundary condition (8). 



III. STATIONARY STATES 

In this chapter we want to explain how the complex 
domain structures found experimentally and from com- 
puter simulations [5,7,13,14] are generated by this simple 
model. 

We denote the electric field profile and the current 
density of stationary states by E* and I* , respectively. 
An easy way to construct the stationary profiles is to 
fix I* , find out the corresponding electric field profile 
{E*}, i = 1, . . . , N , calculate their voltage as a function 
of I* 



The field profile must fulfil the equation [14]: 



1 N 

8 = 1 



(12) 



E*_ 



E* + v 1 



r 



<E*) 



f(E*,n (13) 



The boundary condition implies 
I 



< E *i) 



c+1 



O /(£?,/*) = (14) 



which has three solutions E* for a known fixed value of 
the current on the interval (1 + c) v m < I* < 1 + c. 

In order to understand the properties of the stationary 
profiles we will now investigate the behavior of the set 
{E*} as a function of E\. At first we restrict ourselves 
to c > and monotone increasing field profiles. To con- 
struct {E*}, we have to invert the function f(E,I*) for 
a fixed value of I* . Its derivative is: 



df(E,P) 
dE 



1 



I* dv(E) 
'v(E) 2 dE 



(15) 



With the restriction to monotone increasing field pro- 
files, we can always obtain E* for E* > E m because 
dv/dE > holds. If E{ < E m , Eq. (14) implies 7* < 1 + c 
and we find that f(E, I*) is strictly monotone increasing 
for all E if 



i>max(4±^n 
v ~ \ v(E) 2 dE J 



(16) 



For our function v(E) this yields v < 0.195/(1 + c). In 
this case the function f(E, I*) (I* fixed) is always invert- 
ible. Then we can find a unique field profile parametrized 
by the point E\ . Since we have three possible solu- 
tions E{ of (14) for each given I* /(I + c) £ (v m ,l), 
there are three different voltages <f> for each value of 
the current in this range. The function <f>(I*) is thus 
three- valued, which means that by inverting it we ob- 
tain an N- or an Z-shaped current-voltage characteristic 
as shown in Fig. 1 for v = 0.05 and v = 0.15, respec- 
tively. Both types can be easily understood: When the 
doping density v is low, Eq. (13) shows that E* ps E* 
holds. Thus, the field profile is nearly uniform and the 
current-voltage characteristics follows the t;(£')-curve as 
shown in Fig. 1 for v = 0.05. This is physically obvi- 
ous as there are few charges present inside the sample. 
For larger values of v the values E* may strongly devi- 
ate from E* with increasing i, if E\ is not a fixed point 
of f(E,I*) (which is the case for c = [14]). Let us 
denote by E (1 ~>(I) < E^ 2 \l) < £'( 3 )(7) the three solu- 
tions of v(E) = I for a given I £ (v m , 1) which are the 
fixed points of Eq. (13). If <f> is small and c > 0, E\ 
is on the first branch of v(E) and the values E* tend 
to £ , ( 1 )(7*). When <j> is larger and E* is located on the 
second branch of v(E), the sequence Ti 1 * leaves the neigh- 
borhood of E^(I* ) and then approaches E^ 3 \l*) on the 
third branch of v(E) if v is large enough. This is shown in 
Fig. 2(a). In this case the voltage is basicaly determined 
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by the fixed point E (3 ~>(I*). Since dE (3 ~>(I)/dI > 0, this 
branch of stationary solutions may exhibit a range of 
positive differential conductance leading to the Z-shape. 
This effect is more pronounced for longer superlattices 
with many wells N and also for larger values of c. 

For larger doping v the condition (16) is violated and 
the function f(E, 7) may not be invertible for some cur- 
rent 7. In this case there can be more than one possible 
Ei + i following a given E{. Then the current- voltage char- 
acteristic can no longer be unambiguously parametrized 
by the point E\. In general f(E,T) has 3 different 
branches for a certain interval of 7* , as shown in Fig. 2(b). 
Let us call branch a that having df/dE > for low 
E, branch [3 has df/dE < 0, and branch 7 again has 
df/dE > but for larger E. 

Let us explain how to construct different stationary 
field profiles for a given value of the current I* . We 
shall assume that the profiles are monotone increasing, 
E* +1 > E* , Vi. First of all, E* maybe located on branch 
7 of f(E, I*), and so will be all successive fields E* . This 
profile will have the largest possible voltage for the same 
I* . Secondly, E\ may be located on branch [3, which 
implies that all successive E* of a monotone increasing 
field profile have to be on branch 7. The corresponding 
voltage is smaller than that of the previously described 
branch but larger than those stationary solution branches 
that we analyze next. 

If E* is on branch a, we may have E* , i = 1, . . . , j — 1, 
(j = 2, . . . , N) on branch a, and Ej either on branch [3 or 
on branch 7. We obtain a different branch of stationary 
solutions for each such possibility. Let us denote by ( j, f3) 
or ( j, 7) the solution branch having Ej either on branch 
[3 or 7 respectively, and E* , i = 1, . . . , j — 1 on branch 
a. In Fig. 2(b) a solution (j,/3) is shown by a dashed 
line and a solution (j, 7) by a full line. Clearly j = 1 
corresponds to the possibilities discussed above. Finally, 
we have one solution where all field values are on branch 
a which we denote by (N + 1,7). In order of increasing 
voltages, we have 

<j>(N+l,-y)(I*) < <t>(N,p){I*) < <j>(N,-y)(I*) < 

■ < <?W)(-n < ^(i, 7 )( r )> 

corresponding to 2N + 1 different stationary solution 
branches with the same current I* . They can be ob- 
served in Fig. 1 for v = 1.0 and 7* = 0.8. Notice that 
the branches (j + 1,7) and (j,/3) coalesce at a current 
7* 6 (1, 1 + c) which is roughly independent of v (see 
Fig. 2(c)). The branches (j,/3) and (j, 7) coalesce at a 
lower current 7 C which decreases as v increases (see Fig. 
2(d)). The current- voltage characteristic curve is thus 
connected as shown in Fig. 1 (y = 1.0). 

The field profile of the solution branch (15,7) i s de- 
picted by the crosses in Fig. 3. One can clearly identify 
two regions 1 < 1 ' < j and j < i < N where the electric 
field E* is roughly constant and close to a fixed point with 
v(E*) pa I* . In between there is a transition layer, the 
domain boundary, consisting of only a few wells. These 



type of states we call domain states. A shift of the do- 
main boundary by one well only changes the voltage as 
long as the transition layer does not extend to one of 
the contacts. As the stationary solutions resulting from 
a one- well shift are very similar, the domain branches 
in the current-voltage characteristics look alike, as can 
be seen in Fig. 1 (v = 1.0). The slope of the different 
domain branches may vary with the position of the do- 
main boundary due to the different conductivity of the 
two domains which is related to the shape of the v(E) 
curve. 

Note that the domain states are not very sensitive to 
the exact type of boundary conditions if two conditions 
are fulfilled: 

1. The boundary conditions must allow for the exis- 
tence of a roughly constant field distribution E* ps 

and E* m E^ probably after a short contact 
layer of some wells. 

2. The domain boundary must be located sufficiently 
deep inside the sample, so that it does not collide 
with the contact layer. 

As v decreases, the solution branches become shorter 
and eventually disappear if I c becomes larger than (1 + c), 
which happens if the inequality (16) holds. The station- 
ary domain structures are seen for narrow current inter- 
vals about I* = 1 for intermediate doping as shown in 
Fig. l(v = 0.3). Another complex feature can be found 
here for larger voltages where extra wiggles appear. They 
occur if E* crosses the value 1, yielding an additional 
maximum in I* . 

Except for the fundamental question of stability, we 
have now understood the morphology of the complicated 
current- voltage characteristic curve shown in Fig. 1, its 
changes with doping and its relation with the electric 
field profile. The very same features occur in the more 
complicated model of Prengel et al. [13,17] as shown nu- 
merically in Ref. [25]. 

So far we have restricted ourselves to monotone in- 
creasing field profiles. Therefore we could only find do- 
main states where the high-field domain is located at the 
receiving contact and the domain boundary is an accu- 
mulation layer. Nevertheless, this is not the full story. 
For sufficiently large v other solutions are also possible, 
even for c > 0. A typical such field profile is depicted 
by circles in Fig. 3. The field starts on branch 7 of the 
function f(E, I*) and first increases with the QW index 
towards the third fixed point £'( 3 )(7*). At a certain QW j 
the field jumps down to either branches a or [3, and then 
decreases down towards the first fixed point £ , ( 1 )(7*). Of 
course v has to be large enough for these jumps down 
to be possible. Thus these field profiles have a high-field 
domain located at the injecting contact and the domain 
boundary separating this domain from a low-field domain 
is a depletion layer. Numerical investigation shows that 
these stationary states are stable and that they can be 
reached from many initial conditions. Lemma 2 tells us 
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that the initial field distribution can not be monotone 
increasing in the well index (like the solutions of the con- 
nected branch discussed before), for otherwise the field 
distribution would stay monotone increasing for all times. 
Thus, these different solutions are not connected to the 
branches dicussed before (having a monotone increasing 
field profile) but form many additional isolated curves 
called isolas [26]. A typical isola is shown inside the 
frame in Fig. 1 for v = 1.0, which is also blown up to 
an enlarged scale for the sake of clarity. 

A special situation arises if c = holds. In this case 
the branches of monotone increasing and decreasing field 
profiles become connected and there appears much ad- 
ditional degeneracy leading to an extremely complicated 
structure as observed in Ref. [14]. 



then determines the possible eigenvalues A. For A ^ we 
now introduce the variable Yi = Xii/j and the parame- 
ters 



h = r 



Then we obtain: 



AE1 



\ + a,iYi-i 
Yi. = - — or 



A + a; + bi 
X 



A + &i 



(21) 



A + a; + bi = 
and Yi-i = — A/a; 



(22) 
(23) 



IV. STABILITY OF THE STATIONARY STATES 

Up to now we have only discussed the existence of sta- 
tionary states, but not their stability properties. First 
of all, several stability properties can be established by 
topological arguments [26,27]. If several branches overlap 
at a fixed voltage, each second branch has to be unsta- 
ble by general reasons. For example the middle branch 
(exhibiting positive differential conductivity) of the Z- 
shaped charactistics in Fig. l(v = 0.15) has to be unsta- 
ble. The remaining branches may exhibit further bifurca- 
tions. In order to elucidate this, we will perform a linear 
stability analysis for the states constructed in the previ- 
ous section. By this method we will prove the following 
statements: 

A. For large doping No exceeding approximately 
e(E m - E M )/(ql) ■ v~ m /(v~ M - v~ m ) we find stable 
domain states. 

B. For very small products Nd(N — 1) the almost uni- 
form states are also stable. 

For a medium range of doping in between these two 
limits we find self generated oscillations as reported in 
Ref. [10,11], and further discussed in the Section V. 

In order to perform the linear analysis of stability, we 
set 

Ei(t) = E* +e xt e i (17) 
I(t)=I*+e xt j (18) 

and obtain 

^.j.^,.^,.,,.,, (lw 

This linear equation together with the boundary condi- 
tion eo = ii determines all e; as a function of A. The 
fixed bias condition 

N 

X> = °> (20) 

8 = 1 



Lemma 3 // bi > holds for all i = 1, . . .N the real 
parts of all eigenvalues have to be negative, i.e., the state 
is stable. Furthermore the current-voltage characteristic 
exhibits a positive slope dl* /d<&. 

We prove this Lemma by contradiction. Let us assume 
that Re(X) > holds with A ^ 0. As b x > we directly 
find that Re(Y\) > and \Y\ — 1| < 1. In order to satisfy 
the voltage condition (20) we conclude that there must 
be at least one Yj with Re(Yj) < 0. This implies directly 
that also \Yj — 1| > 1 must hold. But we find: 

, Y k(^-i-i)-M k||>»-i-i| + N 

1 1 ' |A + & 8 - + ai | \bi + ai \ 

(24) 

Given that \Y\ — 1| < 1, the last equation implies that 
\Y - 1| < 1 for all i. Thus, the case Re(X) > 0, A ^ is 
excluded. 

For A = we obtain 
Mi = j and (bi + a^e; = j + . (25) 

Therefore all e; have the same phase as j and the volt- 
age condition (20) cannot be satisfied unless e; = 0, Vi 
and j = 0, which is the trivial case. In conlusion A = 
is not an eigenvalue. Furthermore A = describes the 
infinitesimal change along the curve of stationary states. 
Eq. (25) tells us, that j / ii > 0, Vi. Identifying dl* = j 
and d(f> = ^2 e i we obtain a positive slope of the current- 
voltage characteristic, i.e., dl* j 'd<f> > 0.» 

Therefore we can conclude that the states having all 
their fields Ei in the positive differential mobility region 
are stable, which coincides with our physical intuition. 
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A. Stability for sufficiently large v 

Lemma 3 establishes that stationary field profiles with 
v'(E*) > 0, Vi are linearly stable. These profiles include: 
(i) trivial ones where all the E{ belong to the same branch 
of v(E), and (ii) profiles where the negative differential 
mobility region is crossed in a single jump. This means 
that for a certain value Ej < 1 (1 < j < N — 1) there 
exists E* +1 > E m with f(E* +1 ,P) = E* . As f(E,P) 
is monotone increasing for E > E m , the necessary and 
sufficient condition for the existence of such a value E\ 
is f(E m ,P) < Ej. This condition yields 



+i 



(E m -E*)v(E m ) < v{I*-v m ) 



(26) 



This inequality is first fulfilled for Ej = 1 and I* = 1, 
which are the largest values of the respective quantities 
for the low field domain. This gives: 



v > 



(E m - 1) v„ 

1 - Vm 



(27) 



If this condition is fulfilled, domain states are possible 
which cross the negative differential mobility region in a 
single jump and must be stable therefore. 

Nevertheless there can be stable states even for smaller 
doping v, as Lemma 3 only yields a sufficient and not a 
necessary condition for stability. For our v(E) curve in- 
equality (27) yields v > 0.32. Indeed we do not find 
any self-sustained oscillations for v larger than the value 
v pa 0.27, which is somewhat smaller than our estima- 
tion. Checking other v(E) curves we have always found 
oscillations up to a doping roughly 15 — 40% lower than 
determined by the bound (27). Thus, the bound is not 
only a sufficient condition but also a reliable rough esti- 
mate for doping above which the oscillations disappear. 

Transforming to physical units we obtain a surface 
charge density per well: 



INn > 



v m e (E m — Em) 
(v~m ~ v~ m ) q 



(28) 



which should be a reasonable approximation for the nec- 
essary doping density. For the model equations from 
Prengel [13] we find IN D > 2.4 x 10 n /cm 2 . Actu- 
ally, oscillations are found in the model up to roughly 
INd ~ 10 11 /cm 2 for the regular superlattice and for 
somewhat higher values for a slight amount of disorder 

If we regard domains with the high-field domain lo- 
cated at the injecting contact the same arguments yield 
a bound 



v > 



{ E m ~ 1) 

1 - v m 



(29) 



for the existence of stable domains which is larger by 
a factor l/v m . For our v(E) curve this corresponds to 
v > 0.97. Indeed we have found stable domains with an 



depletion layer for v = 1.0 as depicted in Fig. 3. This in- 
dicates that this type of domain only appears for larger 
doping. This might explain that two different locations 
of the high-field domain have been reported in the lit- 
erature. In Ref. [4] it is found to be located at the in- 
jecting contact for a superlattice with INd = 8.75 x 10 11 
cm -2 while in Ref. [28] the high-field domain is located 
at the receiving contact for a different superlattice with 
IN D = 1.5 x 10 11 cm" 2 . 



B. Stability for sufficiently small v 



Now we want to show that for sufficiently small dop- 
ing, the (connected) branches of stationary solutions are 
stable. Then no self-sustained oscillatory branches bifur- 
cating from them can exist. In order to do this, we note 
that for very low voltages the stationary state is stable 
as indicated by Lemma 3: all field values of this state are 
in the range < E* < 1. We now increase the voltage 
and study whether an instability may occur by checking 
whether it is possible to have A = iu> with u> > for some 
<f>. (The case A = yields the saddle-node bifurcation at 
the point with d<f>* / ' dP = 0, which causes the switching 
to another branch of the Z-shaped characteristic but typ- 
ically does not generate any oscillatory behavior.) In the 
appendix we show that this is possible only if 



(N — l)v > min 



holds, with 



7T (Vj - VCl) 

4ci 



Vl Cl 



2C(I* -vi)c 2 



(30) 



vi = min{i; m , (1 + c) 1 } 

d = I*max El <E<E h \dliLv(E)/dE\ 

c 2 = rmax El <E<E k \d 2 lnv(E)/dE 2 \ 



C 



(N - l>ci 



exp 



(N - l>ci 

vi - VCl 



Here Ei,Eh denote the minimal and maximal values of 
the field for the stationary field profile. Note that for 
small v we find C —> vi/(vi — vci) and furthermore the 
terms vc\ become neglegtible, so that the right side nei- 
ther depends on v nor on N but only on the shape of 
v(E) and the parameter c. 

If v is smaller and the inequality (30) is violated, no 
bifurcating oscillatory branches can issue forth from the 
steady state which is thereby stable. 

The bound (30) is far too small due to the rough esti- 
mations made during its derivation. Therefore the num- 
ber itself should not be used for quantitative investiga- 
tions. Nevertheless we now have shown that the sta- 
tionary states are stable for low doping and that in the 
limit of long superlattices the critical doping decreases as 
1/(N-1). 
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C. Consequences for the continuum limit 

With respect to the continuum limit, N —> oo, v — ► 
0, L := Nv < oo we directly find that there exists a 
minimal length L m such that the stationary state is lin- 
early stable if L < L m . This lower bound is given by the 
Eq. (30) with v = 0, C = 1 and it can be derived directly 
from the equations valid in the continuum limit, as we 
shall report elsewhere [22]. Eq. (30) and similar bounds 
derived for other boundary conditions constitute an ex- 
plicit form of the well-known NdL criterion of the Gunn 
effect [29]: The dimensionless length L (proportional to 
doping times the semiconductor length [30]) has to be 
larger than a certain number for the stationary solution 
to be unstable. 

Obviously the upper bound in the doping v (for the 
absence of the oscillatory regime) does not exist in the 
continuum limit. The discreteness is essential for the 
field distribution to jump from the low-field region to the 
high field region without any fields exhibiting negative 
differential velocity in between, (which stabilizes the field 
distribution). This explains that these stable stationary 
domains can not be found in the usual Gunn diode. 



V. SELF-SUSTAINED OSCILLATIONS IN FINITE 
SUPERLATTICES 

Here we report the results of numerical simulations of 
the self-sustained oscillations of the current in the dis- 
crete model for different values of the dimensionless pa- 
rameters. We have solved Equations (6) and (7), with 
the boundary condition (8) and an appropriate initial 
field profile. 

The simulations yield self-sustained current- 
oscillations for dopings v where the stationary field profile 
becomes inhomogeneous and the middle branch exhibits 
currents that are significantly larger than v(<f>). E.g., 
for the situation from Fig. 1 the oscillations are found 
for dopings v > 0.1. Then the middle branch generates 
oscillatory behavior in a certain range \<\> a ,$u[- We find 
<f> a > 1 and that <j) u is an increasing function of c (due to 
the more pronounced Z-shape as discussed in the third 
section) which can become larger than E m . In the latter 
case bistability between the oscillatory behavior and the 
lower branch of the Z-shaped characteristic, which is sta- 
ble for <f> > <f)p pa E m , occurs if <j) u > <f)p holds. A similar 
type of bistabilty can be found, when the wiggles in the 
static I* (<f>) characteristic appear (Fig. 1, v = 0.3) and 
the upper branch is unstable against oscillatory behav- 
ior while the lower branch is stable. As discussed in the 
last section the oscillations completely vanish for dopings 
where the domain branches are fully developed. 

To clarify the nature of the oscilations Fig. 4(a) shows 
the field profiles at different times of a given period of 
the current oscillations for a 50- well SL. We can iden- 
tify a field profile consisting of two field domains at time 



(1). The domain wall, which is a charge monopole con- 
taining an excess charge of electrons > 1), moves 
towards the receiving contact in time (2-4) where it dis- 
appears. During this process a new domain wall emerges 
slowly as can be seen at the times (3) and (4), so that 
this process is repeated in time. Note that two different 
monopoles are present at the same time in a certain part 
of the period. Between the domain walls regions with 
a more or less homogeneous electric field are observed. 
We define the field of the domain to be the electric field 
E{ at the position where n 8 - has a local minimum, i.e., 
the variation in the field is minimal. In Fig. 4(b) these 
fields of the different domains are depicted by full lines 
as a function of time. Additionally we have shown the 
values of E (1 ~>(I) < E (2 ~>(I) < E (3 ~>(I) (which are the 
three solutions of the equation v(E) = I(t)) as a func- 
tion of time (dashed lines in the figures. Note that the 
solutions E^ and E^ disappear for I > 1). We find 
that the fields of the domains mainly follow the values 
E^\l) and E^ 3 \l). The difference corresponds to the 
displacement current. Monopole recycling and motion is 
clearer for longer SLs as shown in Fig. 5. Notice that 
now the field on the domains follows closely the dashed 
lines, as the displacement current is smaller due to the 
lower frequency. 

The mechanism of the oscillation can be understood 
as follows: Let us start at the time marked by (1) in 
Fig. 5. At this instant we assume that the field profile 
consists of a low field domain with E\ E^\l) next 
to the beginning of the SL and then a high field domain 
with Eh ~ E( 3 \l). Within the low field domain there 
is a tiny inhomogeneity (contact layer) close to the in- 
jecting contact, which is due to the boundary condition. 
The constant voltage condition implies 

N 

El (t) = E h (t) - —(E h (t) - <j>) (31) 

The position j(t) of the domain wall moves to the right 
with a certain velocity v mon /i>, where v mon is always sig- 
nificantly less than 1. Eq. (31) implies that Ei, Eh, and 
I must increase with time so as to fulfil the fixed voltage 
condition. As long as the low field domain is in the region 
of positive differential velocity (Ei < 1), the field profile is 
stable and the contact layer follows adiabatically the cur- 
rent. At to ~ 129 the field E\ becomes larger than 1 and 
the contact layer loses its stability and injects charge into 
the sample. This charge accumulation 6n travels with a 
velocity of the order of v{E\)jv k, 1/v and is amplified 
in time via 

d6n(t) , , „ . „ . . . 

—LL = - v >( El )6n(t) (32) 

with an initial value 6n(to) ~ c. At time (3), to + 7/, 
its size is sufficiently large so that two different field do- 
mains on both sides can be identified and a new charge 
monopole is born. This charge monopole travels than 
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with the velocity v mon jv mentioned above and it sharp- 
ens as it travels; see the point marked by (4) in the fig- 
ure. This stage lasts until the fields before and after the 
monopole reach the values E = E^\l) and E = E^ 3 \l), 
respectively, and we are back in situation (1), having 
completed one period. Mathematically, this behavior can 
be well described by an asymptotic analysis in the con- 
tinuum limit (v — ► 0, Nv = const) [32]. 

Close to the Hopf bifurcation described below we also 
find a different mode of oscillation. We no longer see recy- 
cling and motion of domain walls. Instead, a domain wall 
remains pinned at a given location and the field values 
at the low and high field domains oscillate in antiphase. 
The shape of the current oscillation is almost sinusoidal 
and the maximum current is clearly below 1 (where the 
accumulation layers were injected from the injecting con- 
tact for the mode discussed before). This mode has been 
experimentally observed in [31]. See Fig. 6. A similar os- 
cillation mode has been found in simulations of imperfect 
superlattices [19]. 

Figures 7 show the variation in magnitude and fre- 
quency of the current oscillations as functions of the 
bias for different doping values and number of SL pe- 
riods. In Fig. 7(a), the oscillatory behavior begins at 
<f> = <j) a ?s 1.100 via a supercritical Hopf bifurcation. The 
amplitude of the oscillation increases with bias following 
a square-root law and the frequency is almost constant. 
At <f> = <j) u pa 1.617 the branch of oscillatory solutions 
disappears via a second Hopf bifurcation. For longer SLs 
with N = 200 and the same doping, Fig. 7(e), or for a 
50- well SL with larger doping, 7(c), the end of the os- 
cillatory branch is different: a limit cycle collides with 
the unstable fixed point from the middle branch of the 
Z-shaped current-voltage characteristics and disappears. 
This bifurcation scenario results in a decrease of the fre- 
quency down till zero, when the collision takes place, 
while the amplitude is unchanged in lowest order. In 
all these examples, there is an interval of bistability for 
4> £ {4>f)i4>io)- The bias interval where the oscillatory 
branch exists, (<f> a , 4>io), shrinks as v decreases, and for 
v < v (v k, 0.073 for N = 50, c = 10" 4 ), there is no 
oscillatory solution. 

It is important to note the dependence of the frequency 
with the bias, which could be used to tune the frequency 
of an eventual device [31]. The frequency might increase 
or decrease with bias depending on the sample parame- 
ters N and v, as shown in Figs. 7(b,d,f). 

In order to understand the different types of behavior 
we calculate the period T p of one oscillation. Let us start 
at the time t = to when the the charge accumulation is 
injected at i = 1. One period is completed at t = to + T p 
when the next charge accumulation is injected. At this 
instant, the position of the charge monopole is given by 
Eq. (31) with Ei = 1 and E h « E (3 ~>(1) which yields 



M := j(t + T p ) 



£( 3 )(1)-^ 
£( 3 )(1) - 1 



N. 



(33) 



Thus the charge accumulation has to travel a distance 
M in the time T p . During the monopole formation time 
Tf , the mean velocity of the charge accumulation is v~ x , 
whereas it is equal to v mon jv (with v mon < 1) for the 
rest of the period, T p — Tf . We therefore have 



M =-Tj + ^^(T p -T J ). 



v v 



This gives 



vM ( 1 
T p = - Tf ( 1 



(34) 



(35) 



In order to estimate Tf , we note that is determined by 
the condition 6n(to + Tf)/6n(to) = a where a ~ l/c. 
Then Eq. (32) yields: 



log(a) 



to +?7 



v'(E,(t))dt. 



(36) 



Now we obtain E l (t) from (31) with E h « £" (3) (1), j(t) = 
M + v mon (t — to) /v. Up to the first order in (t — to), this 
yields: 



Ei(t) 



1 Nv mon (E^(l)-^) 



vM 



Linearizing v'(Ei) ps — \v" (l)\(Ei — 1), we obtain 



T 



J 




21og(a) 



,(£(3)(1)-1) 



(37) 



(38) 



Now M decreases as <f> increases, as can be seen from 
Eq. (33). As v mon < 1, Eq. (35) indicates that there are 
two competing mechanisms influencing the dependence 
of T p with <j). In general the monopole formation time Tf 
is negligible compared to T p (<f>) in the following cases: (i) 
for long SLs (M in Eq. (35) is then large), (ii) when the 
values of c are large (i.e., a small), (iii) for large doping 
v. In these cases, T p (<f>) should be a decreasing function. 
We find that the frequency increases with <f> for large val- 
ues of the quantities N, v, c while it decreases otherwise. 
These behaviors are illustrated in Fig. 8. 

Until now we have studied Equations (6) and (7), with 
the boundary condition (8), c, positive. Now we are go- 
ing to describe what happens if — 1 < c < 0, that is, there 
are less electrons in the first well than the doping density 
(m — 1 = c < 0; see Fig. 9). 

As in the case of positive c, the oscillations are due 
to the generation, motion and annihilation of domain 
walls, connecting domains (which are regions of almost 
uniform electric field approximately given by the zeros of 
I — v(E)). The difference is that now the electric field 
profile is monotone decreasing with the QW index: the 
high-field domain is close to the beginning of the SL and 
the low-field domain extends to the end of the SL. The 
domain walls are now charge-depletion layers, having less 
electrons than the doping density. Let us describe one 
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period of the current oscillations for a long SL, such as 
that in Fig. 9 with N = 200. We will assume that initially 
(point marked with 1 in Fig. 9(a)) the field profile has two 
domains connected by a domain wall (E = E^ 3 \lo) to 
the left of the domain wall and E = E^(Io) to the right 
of the domain wall, with an initial value of the current 
£ ( w m;l))- The domain wall is approximately cen- 
tered at j/N = Y = (<f> - (I ))/(E(. 3 \l ) - EW(I )), 
so that Equation (7) holds. The domain wall then moves 
towards the end of the superlattice a with speed close 
to the instantaneous value of the current. The current 
decreases until a certain minimum value slightly smaller 
than v m . Then a new high-field domain is created (close 
to the beginning of the SL) and the current rises sharply 
as the two domain walls move toward the end of the SL. 
When the current is near its maximum value, the old 
domain wall disappears as the values of the field on the 
intermediate and rightmost domains coalesce. Then one 
period of the oscillation is completed. All these features 
can be understood by means of an asymptotic analysis 
to be reported elsewhere [22]. 



stationary domain structures is higher for profiles having 
depletion layers instead of accumulation layers between 
the different domains. 

The transition from oscillatory behavior to a stationary 
state may occur via two different scenarios. For certain 
parameter values, we have found that the stationary and 
oscillatory solution branches coalesce by means of a su- 
percritical Hopf bifurcation. In this case the amplitude 
of the oscillations drops to zero continuously at an almost 
constant frequency. On the other hand, the oscillatory 
solution may collide with the homoclinic orbit of a saddle 
point. In this case the oscillation frequency drops to zero 
while the amplitude remains finite and it does not change 
very much. The latter case has also been found numeri- 
cally for different parameter values and it resembles the 
behavior found experimentally under illumination [10]. 
Depending on the charge at the first QW, the doping 
and the shape of the velocity curve, the frequency of the 
oscillation may depend in different ways on the dc voltage 
bias. This is also corroborated by experiments [10,11,31]. 



VI. CONCLUSIONS 

In this paper we have shown how the complex sta- 
tionary current-voltage characteristic exhibiting domain 
branches is generated continuously as the doping in- 
creases. For low doping the characteristic follows the 
local v(E) relation. If more charges are present, the 
characteristic becomes Z-shaped. When the doping is 
even larger, wiggles appear. For each doping the char- 
acteristic is connected, and the field profiles of all its 
different branches are monotone. The different discon- 
nected branches observed experimentally correspond to 
the stable solution branches of the full stationary current- 
voltage characteristic. It would be very interesting to in- 
vestigate whether it is possible to stabilize the unstable 
branches so that the full characteristic could be observed, 
as in the case of the double-barrier resonant-tunneling 
diode [33,27]. Additionally, for large doping there ex- 
ist isolated branches (isolas) on the full current-voltage 
characteristic having non-monotonic field profiles. 

The stability analysis shows that the almost uniform 
field profile is stable for low doping. The critical dop- 
ing above which time-periodic oscillations of the current 
appear is inversely proportional to the sample length for 
fixed superlattice parameters. This is the same situation 
as in the famous NpL criterion for the Gunn Diode. For 
yet larger doping the time-periodic oscillations of the cur- 
rent disappear: there is an upper critical doping above 
which there appear stable stationary solutions with two 
electric field domains (separated by an abrupt domain 
wall extending almost one period of the superlattice). 
Obviously, this is not possible for the conventional Gunn 
Diode due to the lack of discretization. It is important to 
mention that the upper critical doping needed to stabilize 
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APPENDIX A: PROOF OF STABILITY FOR 
SUFFICIENTLY SMALL v 

Here we prove that A = iu> with u> > can be an eigen- 
value of the linearized system (19) only if v is sufficiently 
large. 

In order to do this, we assume that a given station- 
ary field profile {Ef\ exhibits A = iu> with u> > and 
derive several necessary conditions for this. Restricting 
ourselves to monotone increasing field profiles, E\ must 
either be located on the first or second branch of the 
v(E) curve, as otherwise the second branch is not reached 
which is a necessary condition for the instability accord- 
ing to Lemma 3. Let us now determine the smallest value 
Ei and the largest value Eh the stationary field profile 
{£"?} may take. E\ is given by the value of E\ on the 
first branch for which the current takes on its minimal 
value, I* = 1, considering that the field must eventually 
take values on the second branch of v(E). Eq. (14) yields 
v(Ei) = l/(l+c) which determines the field E\ < 1. Eh is 
given by the largest value that E C N can take on the third 
branch of v(E). Noticing that I* < 1 + c in Eq. (14), we 
can adopt Eh as the solution of v(Eh) = 1 + c from the 
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third branch. Thus, E\ and Eh depend only on the v(E) 
curve and the parameter c but not on the field profile 
{£"?}. For sake of convenience we introduce the follow- 
ing quantities: 



vi 



1 



l + c 



I* max El <E<E h 
c 2 := I* max E ,<E<E h 
Then we have 



dlnv(E) 



dE 
d 2 Inv(E) 



dE 2 



cii > — and \bi\ < c\ Vi 
v 



(Al) 
(A2) 
(A3) 

(A4) 



In the following we will assume that v is so small that 
the function f(E,I*) (I* fixed) is always invertible and 
furthermore 



ai + bi > v i jv - ci > , Vi 



(A5) 



holds. 

For A = iui, ui > 0, we have ReYi > 0, and |Yi - 1| < 1. 
As previously explained in the proof of Lemma 3, there 
must be a Yj such that ReY 8 - > 0, i = 1, . . ., j — 1, and 
ReYj < in order to fulfil the voltage condition. We are 
going to prove the following result: 

Lemma 4 Let j > 1 be the index that satisfies ReYi > 0, 
i = 1, . . . , j — 1, and ReYj < 0. 
(a) If u> < c\, we have 



4A Cl 



where A is the maximum of the expressions 



j 1 
1W ^ a,; + 



(j - k + l)v ^ ^ + bi 

i — k 



(A6) 



(A7) 



for k = 2, . . .j. 

(b) If u> > c\, we have 

(j -\)v> 



VlCl 



where 



B 



2B(I* -v,)c 2 



- 1) (i - 1) 

H := max E ,<E,<E h 



with 



\a,i + bi + ici\ 



(A8) 

(A9) 
(A10) 



Proof: 

(a) Let ijj < c\. 
In order to prove (A6), we consider how the argument 
of the complex quantity Yi is varying with i. 



Yi Y ai y K . ^ ,-r 



(All) 



Therefore we get: 



- fa = arctan ( ) (A12) 



Furthermore we have: 



bi = arctan — 



bi + ai+ujImiY; L ) 



(A13) 



Straightforward calculations starting from eq. (22) yield: 
uj 2 + ciiUjIm(Yi-i) + ai(bi + ai)Re(Yi_i) 



Re(Y) 



ui 2 + (bi + cii)' 1 



Re(Yi) = —^—lm(Yi) + -^—ReiY^) (A14) 
®i t eii Oi -\- cii 

By definition of the index j (j < N), ReYj < and 
ReYj_i > 0. Then these equations indicate that Im Yj < 
and ImYj_i < 0. Thus the transition ReYj_i > —> 
ReYj < occurs across the angle <f> = — 7r/2 as we have 
-7r/2 < <f)j_i < and —w < <f>j < —w/2. 

We introduce the index j' which is defined by the re- 
lations — 7r/4 < <f>j'-i and <f>i < — 7r/4 for i = j' , . . 
Obviously, we have ImYj- < and therefore bi + a; — 
> for i = f,ji + l,...,j. Using (A14) 
and ReYj_i > (for all i < j), we obtain ReYj- > u> 
ImY 8 '/(a 8 ' + bi). This yields for i = j' , . . . , j: 



arctan 



< arctan 



arctan 



|y.i 2 



, . Im(Y,) 

bi + Cii - U |y.| 2 



u> 2 Im(Yj) 

(t,+q,)|y,| 2 



bi + cii -ui : 



, ai + bi 



< 



\Y,\i 
UJ 

ai + bi 



(A15) 



Now we have to distinguish two different cases: 
i) f > 2: By summing the inequality (A15) from i = j' 
to i = j and then taking into account the definitions of 
j and j' , we find: 



7T 7T 

"4 + 2 <^"-i 



< 



E 



ai + bi 
L0(j - j' +l)vAj, 



(A16) 



where definition (A7) has been used. The property 
ui < c\ then implies 
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3 ~ 1 > 3 ~ f + 1 > 



TTVj 



AAjivci 



(A17) Equation (13) for the stationary state now yields < 



ii) j' = 1: This means that <f>i < — 7r/4 and according to 
to Eq. (A13), ijj < —b\ = |&i|. Now we sum the inequal- 
ity (A15) from i = 2 to i = j and then use the expression 
arctan a; > ttx/A for < x < 1, thereby obtaining 



7TU) ( LO \ ( b\ \ 7T 

< arctan = arctan — M 

-46i \-bJ W 2 



< 91 - 9j < 



w v ( j - 1) A 2 



vi 



Therefore we find with Eq. (A4): 

TTV l TTV l 

3 ~ 1 > \ ~ r > 



Av A 2 \bi\ - AA 2 vci 



(A18) 



(A19) 



Putting together (A17) and (A19) we obtain the in- 
equality (A6). 

(b) Let now u> > c\. 
To prove the inequality (A8), we shall define the auxiliary 
functions 



Zi — Y, 



X 



(A20) 



X + bi 

These functions solve the following discrete equation 

(A + a,i + bi) Zi - di Zi_i = Xdi [ — — 1 

VA + 6i_i X + bi 

(A21) 

with the boundary condition Z\ = 0. The solution of 
this problem is 



X (b k - 6 fc _i) 



n 



k=2 (A + (A + b k ) fj- k X + bi + ai 



(A22) 



As all bi are real quantities, we have \iu> + bi\ > u> and 
obtain the following inequality for \Zj \ by using the pre- 
ceding formula with A = iu>, to > 0: 



\z J \<-j2\b k -b k _ 1 \ n 



k = 2 



\bi + ai + iu>\ 



(A23) 



Now we have \Zj \ > - ReZj = - ReYj + to 2 /{to 2 + b 2 ) > 



ijj /(w + bj), where ReYj < and the definition of Z n 
have been used. This inequality together with (A23) 
yield 



< 



u> 2 + b 2 

3 k = 2 



$>*-&*-ii fi 



i = k 



\bi + ai + iu\ 



(A24) 



We now estimate the right side of (A24). The defini- 
tion (21) of bi and the mean value theorem yield 



\b k -b k . 1 \<c 2 \El-El_ l \. 



(A25) 



El - El 



(I*/a k - v), so that < E% - E r k _ 1 < 



J k ^k-l 

v(I* jv\ — 1). Thus we can write 



\b k - <v\ 1 



c 2 ■ 



(A26) 



On the other hand, as we are considering the case u> > c\, 
we find that 

a,; a,; 



< 



\bi + ai + iuj | \bi + ai + ici \ 



(A27) 



according to (A10). Inserting (A26) and (A27) into 
(A24), we obtain 



U) 2 + c{ 



<vc 2 



vc 2 



\k = 2 



I* \ ^~ 

1 n 

vi J H 



j-k+i 



1 



(A28) 



Since u> > c\, u> 3 /(u> 2 + c 2 ) > ci/2. Inserting this into 
(A28), we obtain (A8). Therefore Lemma 4 is proved.* 

Lemma 4 yields necessary conditions for the instability 
of a given stationary field profile {Ef } corresponding to 
a fixed bias. We would like to obtain a general condition 
on v, which should only depend on the v(E) curve and 
the parameter c, but not on the specific stationary field 
profile. This can be achieved by the following considera- 
tions: 



A k 



(3 



i 



l)v ^—j ai + bi 

7 i — k 



< 



Vl 



V} - VCl 



(A29) 



Therefore we have A < «;/(«; — vc\) from Eq. (A7), which 
inserted into the inequality (A6) gives 



(j - l)f > 7T 



4ci 



From the definition (A10) we obtain 

a; 



fi < max El <E,<E 
This yields 



ai - ci 



< 1 



Vl - VCl 



(A30) 



(A31) 



B < 



< 



(j - l)vci 

Vl 



1 



VCl 



exp 



(j - l)^ci 
ted in 

(j -\)v> 



vi - VCl 
(j - l)vci 
vi - VCl 



to be inserted in (A8). The result is 

VlCl 



2C(I* -v,)c 2 



C, (A32) 



(A33) 



We now use the obvious inequality N > j in (A30) and 
(A33) thereby obtaining the condition (30) as a necessary 
condition for oscillatory instability of the steady state. 
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FIG. 1. Current-voltage characteristics for c = 0.01 and 
N = 20 and different values of v. The full line denotes the 
states where the electric field E t is strictly monotone increas- 
ing in i. For v = 1.0 there appear additional branches with 
non-monotonic field profiles E t . They are isolated from the 
stationary branches having monotone increasing field profiles. 
We have shown one such branch, which has also been blown 
up for the sake of clarity. The dotted line is the v(E) curve 
used throughout this paper. 
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FIG. 2. f(E,I) and a trajectory E t indicating decreasing 
i for various doping densities and currents, (a) v = 0.15, 
I = 0.8. (b) v = 1.0, I = 0.8. (c) v = 1.0, 7 = 1.008. (d) 
v = 1.0, 7 = 0.599 = L. 
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FIG. 3. Different stationary electric field profiles for 
c = 0.01, v = f.0, <j> = 1.2, and N = 20. The crosses mark a 
state of the connected branch from Fig. 1(V = 1.0) while the 
circles mark a state belonging to the isolated branch. 
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FIG. 4. (a) Electric field profiles at different times during 
one period of the current oscillation depicted in the inset, (b) 
Time evolution of the electric field values in the left, mid- 
dle (when it exists) and right domain of the SL. The cor- 
responding values of (I(t)), i = 1,2,3, are represented 
with dashed lines. Parameter values are <f> = 1.25, c = 10 -4 , 
v = 0.1, and N = 50. 
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FIG. 5. Same as in Figure 4 for TV = 200. 
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FIG. 6. (a) Same as Fig. 4(a) for N=50 and <f> = 1.61 close FIG. 8. (a) Frequency vs. bias diagrams for different values 

to the Hopf bifurcation, (b) Charge density profiles corre- of the dimensionless doping v and c = 0.001. (b) Frequency 
sponding to the electric field profiles shown in (a). vs. bias diagrams for different values of c and v = 0.1. N = 50. 
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FIG. 7. (a) Stationary current-voltage characteristic (full 
line) with maximum and minimum of the oscillating current 
(dotted lines) N = 50. The oscillatory branch begins at 
4> a ftj 1.100 and ends at <j> w ftj 1.617; the interval of bistabil- 
ity begins at <j>b ~ 1.487. (b) Fundamental frequency of the 
current vs. the average electric field (bias) for N = 50. (c,d) 
Same as (a,b) but now v = 0.15; <j> a ftj 1.052, <j>b ~ 1.556, 
4>u « 1.797. (e,f) Same as (a,b) but N = 200; 4> a « 1.02, 
<f) b « 1.622, <f) w « 1.868. In all cases c = 10" 4 . 
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FIG. 9. (a) Current density versus time when c = —0.01. 
Charge density (b) and electric field (c) profiles during one 
period of the current oscillation. The numbers in (b) and (c) 
correspond to the times marked in (a). Parameter values are 
4, = 1.25, v = 0.35, and N = 200. 
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